#ifndef INITIALIZE
#define INITIALIZE

#include <math.h>
#include "kernel.h"
#include "reduce.h"
#include "framework.h"

__global__ void initializeArrays(float* devData, int devDataPitchInFloats, int nPoints, int nDimension, float* devKernelDiag, float* devAlpha, float* devF, float* devLabels) { 

    //compute index
	int bx = blockIdx.x;
	int tx = threadIdx.x;
	int index = bx * blockDim.x + tx;
	
    //compute diagonal elements of kernel matrix (each thread does a different one)
    //compute its own F (in this case, initialize to the value of the label)
    //set alphas to 0
	if (index < nPoints) {
		devKernelDiag[index] = selfKernel(devData + index, devDataPitchInFloats, devData + (nDimension * devDataPitchInFloats));
		devF[index] = -devLabels[index];
		devAlpha[index] = 0;
	}
}

void launchInitialization(float* devData, int devDataPitchInFloats, int nPoints, int nDimension, float* devKernelDiag, float* devAlpha, float* devF, float* devLabels, dim3 blockConfig, dim3 threadConfig) {
  
    initializeArrays<<<blockConfig, threadConfig>>>(devData, devDataPitchInFloats, nPoints, nDimension, devKernelDiag, devAlpha, devF, devLabels);
  
}

__global__ void takeFirstStep(void* devResult, float* devKernelDiag, float* devData, int devDataPitchInFloats, float* devAlpha, float cost, int nDimension, int iLow, int iHigh) { 
                                 
  float eta = devKernelDiag[iHigh] + devKernelDiag[iLow];
  float* pointerA = devData + iHigh;
  float* pointerB = devData + iLow;
  float* pointerAEnd = devData + IMUL(nDimension, devDataPitchInFloats); 
  float phiAB = kernel(pointerA, devDataPitchInFloats, pointerAEnd, pointerB, devDataPitchInFloats);
	
	eta = eta - 2*phiAB;
	//For the first step, we know alpha1Old == alpha2Old == 0, and we know sign == -1
	//labels[iLow] = -1
	//labels[iHigh] = 1
	//float sign = -1;
 
	//And we know eta > 0
	float alpha2New = 2/eta; //Just boil down the algebra
	if (alpha2New > cost) {
		alpha2New = cost;
	}
	//alpha1New == alpha2New for the first step

	devAlpha[iLow] = alpha2New;
	devAlpha[iHigh] = alpha2New;
	
	*((float*)devResult + 0) = 0.0;
	*((float*)devResult + 1) = 0.0;
	*((float*)devResult + 2) = 1.0;
	*((float*)devResult + 3) = -1.0;
	*((float*)devResult + 6) = alpha2New;
	*((float*)devResult + 7) = alpha2New;

    
    cuPrintf("devResult[0]: %f\n",*((float*)devResult+0));
    cuPrintf("devResult[1]: %f\n",*((float*)devResult+1));
    cuPrintf("devResult[2]: %f\n",*((float*)devResult+2));
    cuPrintf("devResult[3]: %f\n",*((float*)devResult+3));
    cuPrintf("devResult[6]: %f\n",*((float*)devResult+6));
    cuPrintf("devResult[7]: %f\n",*((float*)devResult+7));
    
}

void launchTakeFirstStep(void* devResult, float* devKernelDiag, float* devData, int devDataPitchInFloats, float* devAlpha, float cost, int nDimension, int iLow, int iHigh, dim3 blockConfig, dim3 threadConfig) {
  
    takeFirstStep<<<blockConfig, threadConfig>>>(devResult, devKernelDiag, devData, devDataPitchInFloats, devAlpha, cost, nDimension, iLow, iHigh);
    printf("Test launchTakeFirstStep\n");
}


#endif

